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Viscous heating can play an important role in the dynamics of fluids with a strongly 
temperature-dependent viscosity because of the coupling between the energy and mo- 
mentum equations. The heat generated by viscous friction produces a local increase in 
temperature near the tube walls with a consequent decrease of the viscosity and a strong 
stratification in the viscosity profile which can cause a triggering of instabilities and a 
transition to sec ondary flows. The pr o blem of viscous heating in fluids was investigated 
and reviewed bv lCosta fc Macedoniol l|200.ll) for its important implications in the study 
of magma flows. 

In this paper we present two separate theoretical models: a linear stability analysis and a 
direct numerical simulation (DNS) of a plane channel flow. In particular DNS shows that, 
in certain regimes, viscous heating can trigger and sustain a particular class of secondary 
rotational flows which appear organized in coherent structures similar to roller vortices. 
This phenomenon can play a very important role in the dynamics of magma flows and, 
to our knowledge, it is the first time that it has been investigated by a direct numerical 
simulation. 



1. Introduction 

In this paper we show that the effects of viscous heating can play an important role in 
the channel flow dynamics of fluids with a strongly temperature-dependent viscosity such 
as silicate melts and polymers. In fact, in these fluids, viscous friction generates a local 
increase in temperature near the channel walls with a consequent viscosity decrease and 
often a rise of the flow velocity. This velocity i ncrease may produce a furt her growth of 
the local temperature. As recently described in lCosta fc Macedoniol l|2003t) . above some 
critical values of the parameters of the process, this feedback cannot converge. In this 
case the one-dimensional laminar solution, valid in the limit of an infinitely long channel, 
cannot exist even for low Reynolds numbers. In channels of finite length, viscous heating 
governs the evolution from a Poiseuille regime with a uniform temperature distri bution at 
the inlet, to a plug flow with a hotter boundary layer near the walls downstream l|Pearsonl 
Il977t IOckendonlll97^ . We will show that when the temperature gradients induced by 
viscous heating are relatively large, local instabilities occur and a triggering of secondary 
flows is possible because of v iscosity stratification. 

From previous results (see ICosta fc Macedoniol2003l and references therein) , we know 
that, in steady state conditions for a fully developed Poiseuille or Couette flow, there 
is a critical value of a dimensionless "shear-stress" parameter Q — f3(^) 2 H 4 /(kfio) 



2 



A. Costa and G. Macedonio 



(see below for the symbols used), such that if Q > Q c rit, then the system does not 
admit solution, whereas when Q < Gcrit, the system has two solutions, one of which (the 
solution with greater temp erature) may be unstable. For finite length plane channels, 
ICosta fe Macedonio! l)2003(l have shown that these processes are controlled principally by 
the Peclet number Pe, the Nahme number Na (also called Brinkman number), and the 
non-dimensional flow rate q: 



with p density, c p specific heat, U mean velocity, H half channel thickness, k thermal 
conductivity, po reference viscosity, f3 rheological parameter (see equation Ij2.1|l ). Q flow 
rate per unit length (Q = UH) and dP/dx longitudinal pressure gradient. 
The characteristic length scales involved are the channel dimensions H (thickness) and L 
(length), the mechanical relaxation length L m = UH 2 p/p , and the thermal relaxation 
length L t = UH 2 pep/k. For magma flows, typically L t / L m 3> 1 and the approximation of 
infinitely long channel (from a thermal point of view) is not valid. For finite length chan- 
nels, when viscous heating is important, starting with uniform temperature and parabolic 
velocity profile at the inlet, the flow evolves gradually to a plug- like velocity profile with 
two symmetric peaks in the temperature distribution. The more important viscous dissi- 
pation effects are, the more pronounced the temperature peaks are, the lower the lengt h 
scale for the development of the plug flow is ( Oc kendonlll979ilCosta fc Macedoniol l2003'l. 

Because of the typically low thermal conductivities of liquids such as silicate melts, 
the temperature field shows a strong transversal gradient. Flows with layers of differ- 
ent viscosity were investigated in the past, for their practical i nterest, and it is known 
that they can be unstable depending on their configurati on ijYibl Il967t ICraikl Il969t 
iRenardv fc Josephlll985tlRenardvlll987tlLi fe Renardvlll999]) . In particular, we find that 
when the viscous heating produces a relatively hot less viscous layer near the wall, there 
is the formation of spatially periodic waves and of small vortices near the wall, similar to 
the w aves and vortices wh ich form in core-annular flows of two fluids with high viscosity 
ratio l|Li fc Renardvlll999h . 

In this paper we focus our investigations to the physical regime that typically char- 
acterizes magma flow, with low Reynolds number Re < O(10 2 ), high Peclet number 
Pe ;g> 1, high Prand tl number Pr 3> 1 and low aspect ratio a r = H/L <C I (see e.g. 
IWvlie fc Listell995l) . 

In §|21 we present the governing equations, in §[21 we analyze the linear stability of the 
base flow given by a lubrication approximation, in §0]we describe the numerical scheme 
and the parameters used for the direct numerical simulation (DNS), then we discuss the 
results obtained from DNS and, briefly, few implications for magma flows. 

2. Governing equations 

We consider an incompressible homogeneous fluid with constant density, specific heat 
and thermal conductivity. The fluid viscosity p. is temperature-dependent and, although 
an Arrhenius-type law of viscosity-temperature dependence relationship is more general 
and adequate to describe, for example, the silicate melt viscosities, for simplicity in this 
study we assume the exponential (Nahme's) approximation: 



where T is temperature, (3 a rheological factor and po is the viscosity value at the ref- 
erence temperature Tq. Although a strong viscosity-temperature dependence similar to 
Ij2.f (I , can be responsible for different types of magma instabilities, there have been only 



Pe^pc p UH/k; Na = p U 2 p/k-, q = p Q/(pgH 3 ) 



(1.1) 



p = po exp[-/3(T - T )] 



(2.1) 
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Figure 1. Sketch of the studied system: coordinates and channel dimensions, 
few studies of them (see e.g. IWvlie k Lis"terlll995lll99^) . 

Here we investigate the two-dimensional flow in a plane channel between two parallel 
boundaries of length L separated by a distance 2H (with H/L <1) and we restrict our 
study to a body-force-driven flow (see figure^, although it is not difficult to generalize 
for pressure-driven flow or up-flow conditions for which the driving pressure gradient and 
the gravity act in the opposite direction, as it occurs for example in magma conduits. 
In these hypotheses, the fluid dynamics are described by the following transport equa- 
tions for mass, momentum and energy, respectively: 

V-v = (2.2) 

dv 

p-^ + pv- Vv= -VP + pg + V-T (2.3) 

«f + pv.™ = *V>r + T(J gi (2-4) 

where p is the fluid density, v is the velocity vector, g represents the generic body force, 
P is the pressure, t is the stress tensor, h is the enthalpy per unit mass, T is the tem- 
perature, and k is the thermal conductivity. The term containing the stress tensor t in 
equation l|2.4|l represents the internal heat generated by the viscous dissipation (Einstein 
notation of summation over repeated indeces is used). In this study, for simplicity, the 
latent heat release due to crystallization is not considered, the enthalpy is simply given 
by the product of a constant specific heat times the temperature. Moreover we neglect 
any possible effect due to the buoyancy. Under these assumptions and considering a New- 
tonian relationship between stress tensor and strain-rate (ry = p(dvi/dxj + dvj/dxi)), 
equations l|2.2|l . (|2.3|l and (|2.4|l can be easily expressed in dimensionless form as: 

duj 

= (2.5) 



dui dui 1 „ dp Id _ e / dui duj 
~dt +U] W 3 = ~F~r 9% ~ dtU + E~ed^ f + 

de 80 1 a 96 Nae- e ( du % d Uj 



(2.6) 



(2.7) 



dt ' <K, PedCjdtj Pe 2 \jK, ' d& 

where t = tU*/H is the dimensionless time, (^1,^2) = (x/H,z/H) are the longitudinal 
and transversal dimensionless coordinates, (tti,^) = {v x /U*,v z /U*) represent the di- 
mensionless field velocities (scaled with the characteristic velocity [/*), 9 = [3{T — T ) 
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Name Symbol Definition Value I Symbol Definition Value 



Reynolds number 


Re, 


pU.H/fio 


4.5 


Re 


pUH/Jl 


119.4 


Nahme number 


Na, 


PpoU 2 /k 


14.4 


Na 


ppU 2 /k 


2400 


Froude number 


Fr, 


U 2 ,/(g x H) 


1.5 


Fr 


u 2 /(g*H) 


412 


Peclet number 


Pe» 


pc p U,H/k 


450 


Pe 


pepUH/k 


7400 


Aspect ratio 


Or 


H/L 


3/100 


a r 


H/L 


3/100 



Table 1. Typical dimensionless numbers. The calculated values on the left side are based on the 
mean Poiseuille velocity Up = pg x H 2 /(3/io). The calculated values on the right side are instead 
based the mean velocity U = (J Q v x dz)/H and mean viscosity ~p — (no J exp( — Q))/Hdz. 



the dimensionless temperature, (<?i, g^) = (Sz/lgl, <7z/|g|) indicate the dimensionless body 
force field (from here on- wards we set g z — 0) and p = P/(pU 2 ) is the dimensionless pres- 
sure (Einstein convention of summation over repeated indeces is used). The meaning of 
the usual characteristic dimensionless numbers is reported in Tabled 
Due to the symmetry of the channel and of the boundary conditions, we investigate 
only half of the channel (0 ^ £2 ^ 1). At the walls the boundary conditions are 
given by no-slip velocity and isothermal temperature: Ui = O = at £2 = 1 and by 
112 = dui/d£,2 = dQ/d^2 = at £2 = 0. At the inlet we assume free flow conditions 
and the fluid temperature to be the same as the wall temperature: Gi„ = 0. As initial 
conditions, the velocity and temperature are set equal to zero. 

Considering the geometry of figure ^ an d the isothermal case without viscous heating 
effects, the Navier -Stokes equations of a vis cous liquid driven by a body force g x admit 
a simple solution ijLandau k. Lifschit3ll994ft : 



d 2 v T dP 
Mo^+^ = - = (2.8) 

In this case, the mean velocity is Up — pg x H 2 / (3pLo). 

From this point on-wards, we use starred symbols to indicate the dimensionless number 
based the characteristic velocity Up, i.e. we set £/» = Up, while the un-starred numbers 
are based on the mean velocity U = (L v x dz)/H, i.e. we set U* = U (sec Table H). 
The parameter values used in the DNS and reported in Table are chosen in order to 
perform the computation in a reasonable time, maintaining the system in the regime 
with Re < 0(1O 2 ), Na > 1, Pe > 1, Pr > 1 and H/L <C 1. To fully simulate the 
flow field evolution when viscous heating effects are very important, there is a need to 
solve all the involved length scales of the problem: from the integral length H up to the 
smallest characteristic length-scale. The smallest scales correspond to a thin layer of the 
order of Gz _1 / 2 (ln Na)^ 1 in which the velocity changes from near zero by th e wall to 
near i ts core value (Gz = Pe x H/L indicates the Graetz number) as shown bv IPearsonl 
in the asymptotic limit of very large Na and Gz. 



3. Stability analysis 

The stability of a full y developed steady plane Couette flow was recent ly re-examined 
hvlYueh Wens! (199$, who improved the results previously obtained bv lSukanek et 
(1973). The plane Couette flow shows two different instability modes: one arising in the 
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non-viscous limit, and the other due to the viscosity stratification. As far as the last 
instability mode is concerned, it was demonstrated that the critical Reynolds number, 
above which the flow becomes turb ulent, decreases as th e Nahme number increases, that 
is as the viscous heating increases llYueh fc Wendll99^ . 

Visc ous heating effe c ts on flow stability have been recently investigated experimentally 
bv IWhite fe Mulled (2000), who have shown that above a critical Nahme number an 
instability appears at a Reynolds number one order of magnitude lower than the cor- 
responding Reynolds number predicted for isothermal flow (in these experiments, the 
authors use a temperature-dependent fluid, i.e. glycerin, and a Taylor-Couette device 
which allows the tracking of the vortices by a laser particle tracer). 

When the viscous heating is relevant (Na 3> 1) and the thermal length is much greater 
than the mechanical one, the temperature profile, which is characterized by a narrow peak 
near the channel wall, is dras tically different from the corresponding profile of a thermall y 
steady fully developed flow l|Pearsonlll977t IOckendonlll979t ICosta fc Macedoniol|2003h . 
Assuming slow longitudinal variations of velocity and temperature, we now study the 
linear stability of a thermally developing flow belonging to the important regime with 
Lf /H = Pe ^> 1, LtJ L m ~ Pr ^> 1, Gz 3> 1 that typically characterizes magma flows 



i Wvlie fe Listerlll995|) . In this regime it is legitimate to use a lubrication approximation. 



For the investigation of the linear stability we use the method of small perturbations 
(normal-mode analysis). The base velocity, temperature, viscosity and pressure fields are 
perturbed by two-dimensional, infinitesimal disturbances. Each variable (ui,0,/x,p) is 
given by a steady part plus a small deviation from the steady state: 



where the overbar symbol indicates the steady part, the tilde the perturbation, and v = 
/i//io = e~ e is the dimensionless viscosity. In the (|3.1|1 . the steady part of temperature 
and viscosity depend on the streamwise coordinate £i while the mean flow is assumed 
not to vary appreciably with £i over an instability wavelength. This means that we study 
the thermally developing flow by making the so called quasi-parallcl-flow approximation 
(u2 ~ 0). I.e. one examines the stability of a model flow having the same streamwise 
velocity profile as the real spatially inhomogeneous flow at the selected spatial location. 
Since we treat the stability of those systems in the limit Pe » 1 and Pr ^ 1, with 
the characteristic length L± much greater than the other t ypical mechanical length scales 
l)Pearsorll977tlOckendodll979tlCosta fc Macedonioll2003|) . this assumption is legitimate. 
In this regime it is also legitimate to a ssume th a t the base flow satisfies a system of 
equations similar to that introduced by IPearsonl <jl977|) . At a fixed distance from the 
inlet, we consider the following steady equations: 




3.1. Linear stability 



«i(£i,6,4) = wi (6) +ui(£i,&,t) 



6(6,6,*) = e(6,6) + e(&,&,f) 
K6,£2,4)=*(a,k) + i>(&,£2,4) 



(3.1) 



Jo uid& = 1 




(3.2) 
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with geometry and coordinate system showed in figure ^ As boundary conditions we 
consider V4 = Q = at 6 = ±1 whereas at the inlet (£1 = 0) we assume a parabolic 
velocity profile and an uniform temperature (0 = 0). Equations (|3.2|) were solved by 
a finite-difference method with an implicit scheme for the integration along direction 
6; the pressure gradient was iteratively adjusted at each step in order to satisfy mass 
conservation. 

In the following, we study the linear stability of the base velocity and temperature profiles 
given by (|3.2|) . Since the variations with 6 depend upon the coupling with th e energy 
equat ion through the viscosity, we consider slow temperature variations with 6 l|Pearsonl 
Il977|) . Substituting {ITT} into the equations (|2.5|l . Q2.6|) . (12.71) . subtracting the base flow 
solutions of 1|3.2|) and linearizing, we obtain: 



dui du 2 
96 + 96 







du 2 
dt 



Ml 



~W 
1 dV 

lted& 
du 2 



Ul 



96 
du\ 

96 

dp 



96 

du 2 

96 



d ui d ui 



9% _ 

v d 2 u\ 



96 



96 



dp v 

~96 + Re 
1 dui dv 

R~e~d& 96 + R~e~d& 
2 dV du 2 V ( d 2 u 2 d 2 u 2 
fodb 96 + Re~ \~d& + ~W 



96 2 



(3.3) 



(3.4) 



1 dv dui 
~RedJ'i~d^ (3 ' 5) 




96 \ _ 9 2 9 

U2 W2) ~W 

du 2 £ 9ui 
96 77 96 , 



9 2 e 



+ 



(3.6) 



96 

2VNa^- ( §L 

96 v^6 

Equat ions l|3.3|) , (|3.4|) , (|3.5|) and l|3.6|l are similar to those analyzed bv lPinarbasi fc Liakopouiosl 
( 1995E) who investigated how a variable viscosity affects the stability of the system. In 
this study we account for t he longitudinal variation of the b ase temperature (99/96) 
which was not considered bv lPinarbasi fc Liakopoulol l|l995() and we also introduce new 
terms on the right side of equation H3.6J1 related to the viscous heating. 
In order to eliminate the continuity equation l|3.3|l , we introduce a perturbation stream- 
function ip: 

dip . dip 

Ul = 96 U2 = "96 (3J) 

Moreover we assume that all perturbations have temporal and spatial dependence of the 
form: 

fe,e,i>) = [<M6), 7(6), 0(6), A(6)] e ia ^~^ 

where a is the wavenumber, c is the complex perturbation velocity and <p, /, 
the disturbance amplitudes. 

Substituting equations 1|3.7|) and l|3.8|) into the (|3.3I) . 13.4J) . I|3.5|l and H3.6(l . and eliminat- 
ing the pressure disturbance term by cross differentiation and subtraction, we obtain the 
final stability equations: 



(3.8) 
, A indicate 



iaRe 



217' 



(7J - c)((j>" - a 2 (p) -u </> = V((j) lv - 2a 2 (f>" + a A (f>)+ 
- a 2 cP') + 77"(0" + a 2 <P) + u'(A" + a 2 A) + 2u" A' + u'" A 



(3.9) 



iaPe 



{u-c)6-(t>Q +cp 



,99 
96 



= {9" -a 2 8) + 2vNa 



A 



(3.10) 
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where for simplicity with u we indicate the velocity base flow u\ and the symbol prime 
indicates differentiation with respect to £2- Viscosity perturbation A can be expressed 
in terms of temperature fluctuations by the Taylor expansion of i|2.1 \ . and neglecting 
nonlinear terms: 

A = -99 (3.11) 

obtaining the two final governing stability equations for 4> and 9. Finally, as boundary 
conditions for l|3.9|) and l|3.10fl . we consider: 

= 0, 4>' = 0, 9 = at £ 2 = ±1 (3.12) 

We note that equation 1|3.9|1 reduces to the classic al Orr-Sommerfeld equation when v = 1 
and the equation 1|3.10[) reduces to that used by IPinarbasi fc Liakopoulod l)l995l) when 
both Na = and dO/d^x = 0. 

3.2. Solution method and stability results 

Classical flow stability problems are usually approached in two ways: temporal and spa- 
tial. In the former case, it is assumed that small disturbances evolve in time from some 
initial spatial distribution. In this case, for an arbitrary positive real value of a, the com- 
plex eigenvalue c — cr + icj and the corresponding eigenfunctions (f> and 9 are obtained. 
If cj = Im(c) is negative then the flow is temporally stable, otherwise it is unstable. 
The spatial analysis is focused on the spatial evolution of a time periodic perturbation 
at a fixed position in the flow. This study requires the solution of a nonlinear eigenvalue 
problem in a, which is assumed complex a = cxr + icti with a prescribed real c = cr. 
The disturbances grow for Im(a) < and decay for Im(a) > 0. 

The choice between spati al and temporal study depe nds on the nature of the flow in- 
stability considered (see iHuerre fc Monkewitall990L for a general review ). Moreover 
quasi-parallel flows may contain different region with different stability characteristics. 
In the present paper, a temporal stability analysis of the profiles at a selected set of 
distances from the inlet, has been performed. This analysis is adequate for studying the 
so-called absolute instabilities (i.e. when the perturbation contaminates the entire flow 
both upstream and downstream of the source location) . 

3.2.1. Temporal stability study 

The problem formulated in the S 13. H is solved using a Chebyshev collocation technique, 
expanding the functions (f> and 9 in series of Chebyshev polynomials of order N. The 
2(N+1) coefficients are considered as unknowns and they are evaluated by the collocation 
technique applied at points £2,4 = cos( jy" 3 ) with i = 0, 1, 2, 3...N — 3 and imposing the 
six boundary conditions H3.12J1 at £2 = ±1- This method allows us to define a system of 
2(A^+1) equations in 2(7V+1) unknowns which can be written as a generalized eigenvalue 
problem of the type Ax=cBx. The final system was solved using the LAPACK routine 
ZGGEV. Typically, setting N = 70 and N = 80 permits a satisfactory convergence in 
the computation of the eigenvalues. In order to test the above described compu tational 
implem entation we compared the obtained eigenvalues in the limit O — * 0, with lOrszad 
l)l97llV s results (considering Orszag's definitions, our c is 1.5 times Orszag's c while 
Orszag's Re is 1.5 times our Re). Table shows that eigenvalu es we calculated for 
isothermal limits are very close to those obtained bv lOrszagMl97lh . 
As far as the base flow is concerned, we considered a fixed distance from the inlet £* 
and a given Peclet number Pe. As shown in figure[2]for Pe = 10 7 , as the Nahme number 
increases, velocity distributions deviate from parabolic profile and dimensionless viscosity 
drops near the walls. 
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Mode Number Eigenvalues by Orszag (1971) Our eigenvalues for = 



1 0.23752649 + 0.00373967 i 0.237526311 + 0.00373795 i 

2 0.96463092 - 0.03516728 i 0.964629174 - 0.03516535 i 

3 0.96464251 - 0.03518658 i 0.964643595 - 0.03518749 i 

4 0.27720434 - 0.05089873 i 0.277207006 - 0.05089868 i 

5 0.93631654 - 0.06320150 i 0.936328259 - 0.06320707 i 



Table 2. Least stable e igenvalues c = cr + ici calculated in this work in the isothermal limit 
compared with the lOrszael ljl971f) 's results for Re = 10 4 and a = 1. N = 70 was set. 




-1 -0.5 0.5 1 -1 -0.5 0.5 



z/H z/H 

Figure 2. Base velocity profiles (on the left) and base viscosity profiles (on the right) at £1 = 100 
for Pe = 10 7 and for Na = 0; 1; 10; 100; 1000. Here velocity profiles are normalized with respect 
to the mean velocity U. 

The stability analysis shows that viscous heating in fluids with temperature-dependent 
viscosity is destabilizing. In fact in the cases studied, for a given Pr there is a critical 
Nahme number Na c above which the flow is unstable at any Re, i.e. the critical Reynolds 
number Re c decreases as the Nahme number Na increases. Two clear examples of this are 
shown in figure |3| where, for different values of Na, the imaginary part of the eigenvalue 
c is plotted as a function of the wavenumber a at a distance £* = 100 from the inlet 
and for (Re — 10 2 , Pr — 10 5 ) and (Re = 10 3 ,JY = 10 4 ), respectively. From these 
plots, it is evident that increasing the Nahme number, the imaginary part of the complex 
perturbation velocity tends to increase until becomes positive. For instance, for Pr = 10 
and Re — 10 2 the flow becomes unstable for Na < 10 3 while at Pr = 10 5 and Re = 10 3 
the flow is unstable for Na > 10. The same behaviour was observed with lower Reynolds 
number where the flow becomes unstable at larger Na. 

Beside the investigation of the role of the Nahme and Reynolds numbers on the flow 
stability, a more deepened parametric study of the effects of the other controlling pa- 
rameters of the problem, such as the Peclet number and the distance from the inlet, 
should be performed. In any case, even considering our preliminary results, it is clear 
that viscous heating effects in fluids with temperature-dependent viscosity are impor- 
tant for the determination of flow instablities and, without their inclusion, the critical 
Reynolds number is generally overestimated. 
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Figure 3. Imaginary part of the complex perturbation velocity vs a: Pr = 10 4 , Re = 10 3 (on 
the left) and Pr = 10 5 and Re = 10 2 (on the right). For both cases a distance from the inlet of 
f* = 100 was set. 



4. Numerical simulation 

In this section, we first describe the numerical scheme, then the numerical parameters 
used to solve the equations described above and, finally, we present the results obtained 
and discuss them. 

4.1. Numerical scheme 

To solve equations l]2.2[l . I|2.3I) and (|2.4I) . a fortran code based on the Finite Element 

Meth od (FEM) with the Streamline-Upwind/Petrov-Galerkin (SU/PG) scheme ((Brooks fc Hughes! 
1982ft was used. The enthalpy equation is added in a similar way as suggested by 
Heinrich fc YiJ l|l988l h The solution method is explicit in the velocity and temperature 
and implicit in the pressure, which is computed solving a Poisson equation. 
The domain of interest f2 is partitioned into a number of non intersecting elements f2 e 
with e=l,2,. . . ,n e , where n e is the total number of elements. In contrast with the usual 
Galerkin method, which considers the weighting functions continuous across the element 
boundaries, the SU/PG formulation requires discontinuous weighting functions of the 
form: w = w + s, where w is a continuous weighting function (the Galerkin part) and s is 
the discontinuous streamline upwind part. Both w and s are smooth inside the element. 
The upwinding functions s depend on the local element Reynolds number (momentum 
equations) and the local element Peclet number (enthalpy equation). The SU/PG weight- 
ing residual formulation of the initial-boundary value problem defined by equations H2.2JI 
and H2.3(l can be respectively recasted as: 



(4.1) 



where s^. is a weighting function which is chosen to be constant within each element, and 
discontinuous across the element boundaries, and 



dvi 



PfJi 



8V; 

a \ dt 

[ ~u ( 9vi dvi 



dn- 

da. 



dwk 



dn- 



dt dxj dxj 



P9-> 



dQ, = I aoi u>k dT 



where is the upwinding function for the momentum equation, <7j 



-PSi 



(4.2) 



(Sij is 



the Kronecker symbol). With we generally indicate the boundary surface where the 
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variable (f> is prescribed. 

In the same way, the weak form of the enthalpy equation l|2.4(l may be written as: 



f f dh dh\ m f /, dT\ dw k m 

Jn V at ox i J Jn \ ox i / ox i 

dh dh_ _ _d_f k dT_\ % dvi_ 

_ dt pVl dxi dxi \ dxi ) u dxj 



E 



(4.3) 
dfl = I qodT 



in 

where, is the upwinding function of the enthalpy equation and go is the heat flux 
through the boundary surface I\. In the present work, the velocity field v and the 
enthalpy fields are linearly interpolated with multi-linear iso-parametric interpolation 
functions using rectangular elements. The pressure field P, instead, is assumed to be 
constant within each element and discontinuous across the element boundaries. 
Equations (|4.2(1 and 14.3J1 yield two algebraic equations which may be combined in the 
following one: 

Ma + Cv + N(v) - GP = F (4.4) 
whilst the continuity equation (|4.1J) yields to: 

G T v = D (4.5) 

In the above equations, the vector v represents the nodal values of the velocity Vi and the 
temperature T, whereas the vector a represents the nodal values of the time derivatives 
of the velocity Vi and of the temperature T, P indicates the pressure field, M is the 
consistent generalised mass matrix, C and N(v) account for the diffusive and the non- 
linear convective terms respectively, F is a generalised force vector, D accounts for the 
prescribed velocity at the boundaries, G is the gradient operator, and G T its transpose. 
Equations l|4. 4J1 and (|4.5(l are integrated in time starting from the velocity and pressure 
fields at t = 0. 

Convergence of the algorithm is assured when the element Courant number Cr satisfies 
particular conditions depending on the element Reynolds and Peclet numbers. Typically, 
to assure the convergence, the e lement Courant number should satisfy the most restrictive 
among the following relations l|Brooks k. Hugheslll982|) : 

Cr ^ 0.8, if 7 = (p e e el ^ 100, 

Cr <mirx(l,7), if 7 = (p^ el < 100 

where Cr = vAt/ Axi, with v, At and Axt element velocity, computational time step and 
computational grid size respectively, and 7 = Re e i or Pe e i represents both the element 
Reynolds and Peclet numbers. Unfortunately, the above convergence criteria can be very 
restrictive, forcing the choice of a very small time step to guarantee the convergence of 
the algorithm. 

4.2. Numerical parameters 

Since viscous friction is greater near the walls (higher gradients), it is convenient to use 
a computational grid finer near the boundaries and coarser towards the centre of the 
channel. 

The computational grid was formed by an uniform horizontal grid size Ax/ H = 8.3 TO -2 
while the vertical mesh size A£ = Az/H is not uniform and consists of three different 
grid sizes Ad, A^2, and A^. The finest grid size A£i = 1.67- 10~ 2 was set near the wall 
where the fields change rapidly and the viscosity is lower, A£ 2 = 3 • 10~ 2 was set in the 




x/H 

Figure 4. Zoom of the bottom-left corner of the computational grid used for the simulations. 
The entire computational domain was discretized with 401 x 38 rectangular elements. The lon- 
gitudinal grid size was chosen uniform Ax/H = 8.3 ■ 10 -2 while three different sizes were used 
to form the transverse grid: ni = 23 elements of size Az\/H = 1.67 • 10 -2 near the wall, ni = 5 
elements of size Azi/H = 3 ■ 10~ 2 in the intermediate region, and 723 = 10 elements of size 
AzijH — 4.67 • 10 -2 in the central part of the channel. 



intermediate region and finally, A£ 3 = 4.67 • 10~ 2 in the central part of the channel. The 
computational grid used for the simulations is shown in figure 

A time step of At = 5 • 10~ 4 was chosen to perform the simulation, with one predictor 
and one corrector iterations per time step, whereas spatial integration was performed 
using a 2x2 Gauss quadrature. All runs were performed in double-precision arithmetics 
on a HP-J5600 workstation. 

From a practical point of view, an estimation of the minimum grid size required to well 
resolve the rapidly varying fields was assumed to be equal to the smallest scale involved 
in the problem <|PearsorJl977|) : A^ 0) = Gz" 1 / 2 ^ Na)" 1 (using some Na and Gz initial 
estimations) , it was then further reduced to guarantee the numerical convergence and in 
such a way that the numerical solution does not appreciably depend on the computational 
grid size. The final dimensionlcss numbers used in the DNS are reported in Tabled 

4.3. Results and discussion 
In this section we describe results obtained by the direct numerical simulation (using 
the FEM code described in § 14.1(1 for the values of the parameters reported in § 14.21 and 
Table^for a channel flow of length 100/3 iJ-unit. 

In this study we relied on the small numerical round-off errors present in any numerical 
simulation to trigger natural modes, but further cases with e.g. selected harmonic or 
random input disturbance should be investigated. 

From these simulations we can see that as the time increases the temperature starts to 
rise in the region near the outlet because of viscous heating. At t = t c w 40 an instability 
is triggered in this region where the dimensionless temperature locally becomes greater 
than w 5 (see figurelSJl. For t > t c , as viscous heating effects become more important even 
in the more internal region, secondary flows appear to organize themselves into "coherent 
structures" as rotational flows. This kind of secondary flow looks like roller vortices which 
seem to move from the region near the outlet towards the inlet (see figure 01 ■ Actually 
this happens because the viscous heating becomes relevant even in the internal region 
and the entire flow becomes unstable. 

Figure show the temporal profile evolution at a given distance (for example at £1 = 
22 that is about 2/3 of the channel length). We can see that O, starting with a flat 
distribution, gradually increases near the wall forming a profile with a maximum at a 
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short distance from the boundary. As time increases, this peak becomes more pronounced 
(8 m0 i < 6) filling, at the steady state, a narrow shell of values at a shorter distance from 
the wall (see figure[7}. As a consequence the dimensionless viscosity profile e , strongly 
decreases in correspondence of the temperature peak, reaching values much samller than 
its initial ones (see figure [TJ. The layer where the viscosity is very small is immediately 
close to the colder layer adjacent to the wall and it corresponds to the region where the 
vortices appear (see figures El . 

The longitudinal velocity profile (scaled by Up), starting with a parabolic distribution, 
evolves toward a plug profile filling, at the steady state, a narrow region of values with 
a plug velocity < 18 (see figure [SJ. Figure E] also shows the evolution of the transversal 
dimensionless velocity profiles which, because of the vortical motions, near the wall, tend 
to fill an onion shape region with the largest fluctuations in corrispondence of the peak 
in the temperature profile. For comparison, the profiles computed using the lubrication 
approximation i|3.2f> are reported in figures and El 

Some of these results could be expected on physical basis. In fact, as a first approxima- 
tion, because of viscous dissipation effects, the flow can be viewed as a two-layer-flow of 
two different viscosity fluids with the less viscous one flowing near the wall. The simula- 
tions perfomcd confirm this limit, in fact when viscous heating form a consistent layer of 
less viscous liquid near the wall, the behaviour of this flow tends to be similar to that of 
a two-layer flow with the more viscous fluid in the central part and the less viscous fluid 
near the wall. This arrangement is common in transporting heavy viscous oils which ar e 
lubricated using a sheath of lubricating water l|.ToseDh et aZ.lll997t iLi fc Rena rdv 1999). 
Experiments and simulations of this two-layer flow type of fluids with high viscosity ra- 
tio predict spatially periodic waves called bamboo waves because of their shape, and the 
formation of vortice s in the region near the wall distributed in the trough of the waves 
1 Joseph et al\\\mi\). 

In our simulations, these features can be seen from figure Eland from figure El where a 
zoom of the flow fields near the channel wall is shown. In fact following the flow isolines, 
a spatially periodic wave can be easily discerned and relatively large vortices, settled in 
the middle of the wave troughs, are also evident. 

More over, similarly to the core-annular flows with high viscosity ratio l)Li fc Renardvl 
1999|, the formation of a mixed profile (with a counter-flow zone) near the wall, leads to 
the appearance of vortices (figure EJl. 

Finally, using the dimensionless numbers reported in Tabled we perfomed a linear stabil- 
ity analysis of the base profiles given by the lubrication approximation (|3.2|) at a distance 
2/3 of the tube length (£f ~ 22) from the inlet. These analysis indicate that the base 
flow is already unstable for Na < 120 even at Re = 120 and, as it is shown in figure El 
the most dangerous mode for this flow has wave number a » 7, corresponding to a wave 
length A ~ 1.1 (in 7J-unit) which appears in agreement with that given by DNS. In fact, 
as it is shown in figure El at the distance Q « 22 from the inlet, a wave length of A ss 1.2 
can be estimated. 

In order to obtain a closer comparison of the linear stability theory with the nonlinear 
results obtained by the numerical code above presented, we simulated the case of a chan- 
nel flow like that previously described, with conditions at the inlet given by the solutions 
of the equations for Pe = 7400, Re = 119.4 and Na = 200 at $1 = 21, free flow 
conditions in the remaining boundary and no-slip conditions at the walls; as initial field 
inside the tube we imposed u.- L = = 0. The linear theory predicted for these profiles 
a growth rate a = aci ~ 0.42 (in U/H-unit). To compare the linear with the nonlin- 
ear regimes, the evolution of the maximum amplitude A(t) with time (in H/U-unit) is 
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Figure 5. Evolution of the dimensionless flow fields. The figures show the simulated streamlines 
with the temperature field as background at different times t. The values reported along vertical 
and horizontal axes indicate the dimensionless distances from the channel wall and from the 
inlet, respectively. The blue colour indicates the lowest dimensionless temperature (G = 0) and 
the dark orange corresponds to the highest temperature (O = 6). The symbol t on the upper 
right corner indicates the dimensionless time t. For visualization reasons the horizontal axis is 
contracted with respect to the vertical. 



plotted in figure I1UI The maximum amplitude growth shown in figure 1101 is given by 
the evolution of streamlines around a short distance from the inlet (0 ^ £i ^ 1.3). As 
initial amplitude we considered the value of A(t) at the time at which small perturba- 
tions due to round-off errors begin to grow. Figure ITUl shows that the initial evolution of 
the perturbation is close to that predicted by the linear theory, and then rapidly starts 
to deviate as amplitude increases. After only less than about one time H/U the linear 
growth completely fails in the prediction of the amplitude evolution. 

Th ese and other preliminary results of the investigations on viscous heating effects 
(e.g. ICosta fc Macedonioll2003^ may help in the understanding of some common phe- 
nomena that may occur during lava and magma flows. For instance, effects of viscous 
dissipation can efficiently enhance thermal and mechanical wall erosion, and can help 
to understand the reasons of the inadequacy of simple conductive cooling models com- 
monly used to describe lava and magma flows. Moreover in volcanic conduits viscous 
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Figure 6. Visualization of the flow structures near the channel wall for t — 195 at a dimen- 
sionless distance from inlet of about 2/3 of the channel length. The blue colour indicates the 
lowest dimensionless temperature (O = 0) and the darkest orange corresponds to the highest 
temperature (0 = 6). 




Theta V 

Figure 7. Temporal evolution from t = (blue colour) to t = 250 (red colour) of dimensionless 
temperature profile O (on the left) and dimensionless viscosity profile v (on the right), at a 
dimensionless distance from inlet of about 2/3 of the tube length (£* = 22). The coding colour 
maps the time evolution from the blue up to red. For comparison, temperature and viscosity 
profiles computed using the lubrication approximation (equations ^3.2p ) are also reported (full 
lines with crosses). 



heating could play an important role on the dynamics of both effusive and explosive 
eruptions, influencing directly or indirectly magma gas exsolution and fragmentation 
ijVedeneeva et aDl2005|) . Since magma flowing in conduits and channels is much hotter 
than the wall rock, another dimensionless number B — /3(Tj„ — T w ), that compares the 
imposed difference of temperature with [3 should be considered (Tj n and T w represent 
inlet and wall temperatures respectively). However, previous preliminary studies indicate 
that in m agma flows, viscous dissipa tion effects can overcomes the thermal cooling from 
the walls ijCosta fc Macedoniol2003^ . Moreover, although by increasing B the peak in the 
temperature profile moves towards the centre of the channel, because of the low magma 
thermal conductivities, the flow behaviour i s not much different from the case with B = 
MacedouirJlifini ISchneiderljl 97r1) . 
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Figure 8. Temporal evolution from t = (blue colour) to t — 250 (red colour) of the profile of 
the dimensionless longitudinal velocity u — ui = v x /Up (on the left) and of the dimensionless 
transversal velocity profile v = «2 = v z /Up (on the right), at a dimensionless distance from 
inlet of about 2/3 of the tube lenght (££ = 22). The coding colour maps the time evolution 
from the blue up to red. For comparison, the velocity profile computed using the lubrication 
approximation (equations 13.21 ') is also reported (full line with crosses). 
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Figure 9. Imaginary part of the complex perturbation velocity vs a: Pe = 7400, Re = 119.35 
and Pr = 62 at a dimensionless distance from inlet of about 2/3 of the channel length (£* = 21). 



4.4. Validity and limits of the model 

We have seen that when viscous heating is relevant, a special class of secondary flows can 
develop in fluids with temperature dependent viscosities even at low Reynolds numbers. 
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Figure 10. Maximum amplitude A(t) versus time (in H /U-unit) on a ln-linear scale. Solid 
line represents theoretical linear growth and crosses represent simulation. 



This kind of vortical structures is locally confined near the walls where there is a large 
viscosity gradient and the viscosity is lower. 

The results obtained are valid in the limit of a 2D model based on the full solution 
of the Navier-Stokes equations although turbulence is generally three-dimensional even 
starting with two-dimensional initial conditions. On the other hand, it is known that 
the growth of three-dimens i onal instabilities may be su ppressed by a strong anisotropy 
ijSommeria fc Moreaulll982l iMessade k k, Moreaull200 ll). Th is anisotropy can be due to 
the presence of a mag netic field llSommeria fe More aulll982i). a str ong rotation and/or a 
density stratification <|Lillvlll972UHopfingerlll987t Ivan Heiistfll993|) . 

As in the isothermal case, the Squire's theorem suggests that for the linear stability 
analysis it is sufficient to consider two dimensional disturbances. Although the case of 
core-annular flows with high viscosity ratio suggests t hat a 2D model is able to describe 
well t he flow features observed during the experiments ((Joseph et "^]ll997tfjTfc Renardvl 
Il999l) . because of the complexity of these non-isothermal flows, the effects of 3D distur- 
bances on a quasi-2D flow should be also investigated in order to understand whether the 
evolution of three dimensional motions could be able to obscure the vortical structures 
described above. In our case, we suppose that the strong viscosity stratification induced 
by viscous heating could inhibit 3D motions and the 2D model we used should be able 
to account for the essential physical properties of the real systems; however only an ex- 
tended 3D simulation can completely confirm this. 

Finally, we note that the numerical scheme we used has a first-order upwind and it needs 
very restrictive conditions and a large computational time in order to be accurate. A 
more efficient scheme should be used to permit a more complete parametric study. 
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Conclusions 

The thermo-fluid-dynamics of a fluid with strongly temperature-dependent viscosity 
in a regime with low Reynolds numbers, high Peclet and high Nahme numbers were 
investigated by direct numerical simulation (DNS) and the linear stability equations of 
the steady thermally developing base flow was studied. 

Our results show that viscous heating can drastically change the flow features and fluid 
properties. The temperature rise due to the viscous heating and the strong coupling be- 
tween viscosity and temperature can trigger an instability in the velocity field, which 
cannot be predicted by simple isothermal Newtonian models. 

Assuming steady thermally developing flow profiles we performed a linear stabilty anal- 
ysis showing the important destabilizing effects of viscous heating. 

By using DNS, we showed as viscous heating can be responsible for triggering and sustain- 
ing a particular class of secondary rotational flows which appear organized in coherent 
structures similar to roller vortices. 

We wish our preliminary results can stimulate further more accurate studies on this in- 
triguing topic, contributing to a more quantitative comprehension of this problem which 
has many practical implications such as in the thcrmo-dynamics of magma flows in con- 
duits and lava flows in channels. 

We would like to acknowledge the anonymous referees who strongly improved the 
quality of the paper with their useful comments. We also thank S. Mandica for his 
corrections and suggestions. 
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